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Abstract  A  Normal  Mode  Analysis  of  the  Chesapeake 
Bay  has  reached  the  phase  of  performing  time-series  anal¬ 
ysis.  Prior  attempts  using  generalized  eigenfunctions  ob¬ 
tained  from  COMSOL  have  led  to  extraction  of  power 
spectra  using  a  partial  domain  extraction  and  a  wavelet¬ 
like  correction  to  the  result.  A  new  approach  has  been  de¬ 
veloped  which  effectively  solves  the  spatial  equivalent  of 
the  Initial  Value  Problem.  By  taking  data  sets  at  a  limited 
number  of  points  (  about  10)  over  extended  periods  of  time, 
the  Chesapeake  Bays  flow  vector  field  can  be  extracted 
from  the  power  spectra.  This  method  combines  the  spa¬ 
tial  modes  derived  from  COMSOL  with  the  generalized 
basis  set  for  time.  By  overlaying  windows  of  the  spectra 
obtained,  a  progression  of  the  spectral  components  is  ob¬ 
tained  for  the  lowest  100  eigenmodes.  This  progression 
can  exhibit  smooth  behavior  which  allows  for  prediction 
of  the  modal  behavior  to  take  place,  effectively  predicting 
the  flow  across  the  entire  Bay  while  only  sampling  the  Bay 
at  10  locations. 
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1  Introduction 

The  study  of  vector  flow  fields  for  fluids  has  a  rich  his¬ 
tory.  A  study  of  the  Chesapeake  Bays  flow  has  been  un¬ 
derway  for  several  years  by  many  different  parties.  Tra¬ 
ditional  approaches  involve  modeling  the  Bay  in  terms  of 
fluid  models  coupled  with  geophysical  effects  such  as  the 
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Coriolis  force.  The  approach  taken  by  the  authors  is  not 
to  model  the  Chesapeake  Bays  flow,  merely  to  character¬ 
ize  its  signal  properly  so  that  a  model  may  emerge  from 
the  data.  By  changing  the  emphasis  from  model  validation 
to  data-driven  signal  characterization,  a  more  focused  ef¬ 
fort  is  being  made  to  understand  the  specific  problems  the 
data  presents,  irrespective  of  what  physical  models  may 
or  may  not  be  relevant.  The  method  of  choice  has  been  to 
calculate  Normal  Modes  using  COMSOL.  For  a  historical 
review  of  Normal  Mode  Analysis  (NMA),  please  refer  to 
the  work  based  on  [Eremeev  et  al.  1992]  [11  [2]  and  [Lip- 
phardt  et  al.  2000]  [3]. 

As  a  brief  reminder,  the  two  studies  leading  up  to  the 
Chesapeake  Bay  analysis  were  presented  at  COMSOL  "05 
and  the  World  Congress  on  Computing  CSC’06  confer¬ 
ences  showing  the  Dirichlet  [6],  [9]  and  Neumann  [7],  [9] 
solutions  computed  via  COMSOL  as  well  as  a  finite  dif¬ 
ference  in-house  code.  The  numerical  aspects  of  this  paper 
were  motivated  by  a  method  for  completing  surface  cur¬ 
rent  velocity  fields  called  Normal  Mode  Analysis  (NMA) 
[Eremeev  et  al.  1992] [1]  [2],[Lipphardt  et  al.  2000][3]. 

In  each  of  the  former  cases  studied  using  NMA  the  end 
goal  was  to  be  able  to  use  Normal  Modes  to“fill  in”  vector 
fields  where  data  is  not  present  and  to  extract  power  spec¬ 
tra  in  order  to  search  for  time-dependent  features.  When 
strong  time-dependent  features  are  present,  short-term  pre¬ 
diction  may  be  employed.  Eremeev  et  al.  studied  the  Black 
Sea,  collecting  data  from  autonomous  drifting  buoys  ( ADB) 
Lipphardt  et  al.  used  this  same  method  to  fill  in  gaps  of  ve¬ 
locity  fields  for  Monterey  Bay  using  HF  radar  for  the  data 
set.  The  Black  Sea  is  a  closed  water  body  roughly  in  the 
shape  of  a  kidney.  Monterey  Bay  is  open  to  the  Pacific 
Ocean  and  has  a  shore  roughly  hemispherical.  Lipphardt’s 
group  extended  the  NMA  by  adding  a  mode  to  account 
for  the  flow  between  Monterey  Bay  and  the  Pacific.  The 
Chesapeake  Bay  has  1 1,684  miles  of  shoreline  but  is  only 
189  miles  long  by  30  miles  wide,  giving  it  a  jagged  shore, 
almost  fractal,  compounded  by  a  large  opening  to  the  At¬ 
lantic  Ocean  at  its  southern  end.  The  length  of  the  Chesa¬ 
peake  Bay  allows  for  both  fresh  water  as  well  as  salt  water 
to  exist  within  its  boundary.  As  a  result,  the  Chesapeake 
Bay  represents  a  difficult  system  to  model.  Figures  1,  2 
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Fig.  1  The  Black  Sea  compliments  of  NASA  World  Wind. 


Fig.  2  Monterey  Bay  from  NASA  World  Wind.  The  open  boundary 
with  the  Pacific  clearly  visible. 


and  5  illustrate  the  varying  levels  of  complexity  that  have 
been  solved  using  NMA. 

The  basic  unit  of  calculation  used  throughout  this  pa¬ 
per  is  the  normal  mode.  Like  the  modes  of  a  guitar  string 
or  an  organ  pipe,  systems  obeying  the  Helmholtz  equation 
and  Dirichlet  or  Neumann  boundary  conditions  will  res¬ 
onate  in  states  referred  as  "normal  modes”.  For  the  Chesa¬ 
peake  Bay,  the  modes  calculated  are  energy  potentials  whose 
gradients  and  curls  of  gradients  correspond  to  the  vector 
current  fields  found  in  fluid  mechanics  (it). 


Fig.  3  Chesapeake  Bay  from  NASA  World  Wind.  The  Atlantic 
ocean  open  at  the  southern  end,  allowing  in  salt  water.  The  north 
end  dominated  by  fresh  water. 

Briefly,  the  formulation  leading  to  the  calculation  of 
fluid  flow  stems  from  the  realization  that  the  vector  fields 
can  be  derived  from  two  scalar  fields,  which  are  the  solu¬ 
tions  to  the  Helmhotz  equation  under  Dirichlet  and  Neu¬ 
mann  boundary  conditions  [4]. 

it  —  V  x  [(w'F)  +  V  x  (nO)].  (1) 

Here  ‘F  is  the  stream  potential  where, 

,  .  ,  t-d'V  <?*F\ 

=  (2) 

and  0  is  the  velocity  potential  where, 

->  (dO  dO\ 

=  =  (3) 

with  (m,v)  representing  the  surface  current  velocities  in 
the  x  and  y  directions  respectively.  The  total  velocity  field 
is  composed: 
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Tt(x,t)  =  X  [a(t)n ltDn  +b(t)n TfjyJ  +  it Source-  (4) 

72=0 

The  non-conserving  flow  of  mass  across  the  boundary, 
it  source,  is  model  or  experimental  data  of  the  fluid  veloc¬ 
ity  at  the  boundary.  The  source  term  requires  input  from 
an  external  reference,  in  this  case,  the  work  of  Tom  Gross 
at  NO  A  A  was  employed  where  source  currents  were  pro¬ 
vided  from  the  QUODDY  finite  element  model  of  the  Bay. 
Although  QUODDY  is  a  model  of  the  Chesapeake  Bay’s 
currents,  it  draws  data  from  multiple  sets  collected  near 
the  mouth  of  the  Bay  (the  Atlantic  ocean  interface)  as  well 
as  selected  points  from  the  interior  of  the  Bay’s  geometry. 
For  a  recent  review  of  efforts  by  NOAA  and  Tom  Gross, 
see  reference  [5], 


2  Time  Series  Analysis 

The  Galerkin  method  is  used  to  calculate  the  amplitude 
for  a  given  normal  mode.  A  key  point  to  this  analysis  is 
that  the  technique  requires  knowledge  of  the  data  over  the 
full  domain.  The  source  term  will  be  subtracted  from  the 
outset,  eliminating  it  from  this  analysis.  Only  focusing  on 
the  x-component  of  the  vector  flow  field,  u(x,t),  and  as¬ 
suming  a  data  set  exists  which  spans  the  entire  x-domain. 


that  the  entire  vector  flow  field  can  ever  be  obtained  at  any 
one  time.  Satellite  data  may  be  able  to  provide  information 
related  to  vector  flow,  however,  in  order  to  extract  use¬ 
ful  time  series  would  require  long-term  surveillance  of  the 
area  in  question.  Satellite  data  is  simply  too  expensive.  At 
the  COMSOL  2006  conference,  a  partial  domain  Galerkin 
method  was  investigated  and  found  to  yield  wavelet-like 
responses  in  the  spectral  domain  if  the  data  did  not  cover 
100%  of  the  domain.  This  effect  is  partially  correctable, 
but  unwanted.  In  the  worst  case  scenario,  the  results  are 
unusable  [8],  Figure  2  shows  the  response  of  the  power 
spectrum  for  the  15th  mode  when  only  50%  of  the  domain 
has  data.  For  a  delta  function,  the  response  is  similar  to  the 
function  sinc(k  x)  —  sm^x> ;  prevelant  in  wavelet  theory. 


a  for  mode  m=1 5  out  of  50  modes 


mode  number,  m 


“(*>*)  =  X 

77=0 

a(t)m  =  j  u(x,t)datau{x)DmdQ 

a(t)m  =  f  X  [a(t)AX)o,n  +Kt)AX)N,rI]  “(4,.^ 

J  71=0 

a(t)m  =  X 

72=0 

n  —  (On  &nm  ^(O«0] 

72=0 

a(0m  =  Snma(t)n . 

a(t)n  =  j  u{x,t)dalau{x)Dndn 

similarly,  for  the  Neumann  coefficients,  bn,  coefficients 
are  calculated  by  integrating  with  the  Neumann  solutions, 
u(x)Nn  over  the  domain. 

Central  to  the  theme  of  this  derivation  is  the  applica¬ 
tion  of  orthonormality  of  the  basis  sets.  The  common  rea¬ 
son  for  not  applying  this  technique  is  that  full  knowledge 
of  the  vector  field  is  required  over  the  entire  domain  in  or¬ 
der  to  obtain  the  coefficients  a(t)„  and  bit),,.  For  a  water¬ 
way  like  the  Chesapeake  Bay,  it  is  impractical  to  assume 


a(t)  ®  u.,  u.,  dQ+b(t)  ®  u.,  u.,  d£2 

\  /n  T  D,n  D,m  \  Jn  T  N,n  D,m 


3  Normal  Mode  Analysis 


Taking  a  step  back,  in  solving  a  ID  wave-like  problem 
with  both  spatial  and  temporal  boundary  conditions  may 
reveal  a  different  approach.  Consider  a  guitar  string  of 
length,  L.  Held  at  both  ends,  clearly  a  Dirichlet  problem. 
Using  the  wave  equation: 

c2V2F(x,t)  =  ^-1F(x,t),  F(x,t)=f(x)  g(t).  (5) 

dt 

Rearrangement  of  the  equation  leads  to  a  separation  of 
variables: 

£2 /MS  W  =  ^2  /(*)*(') '  (6) 

fxx  _  Sit  _  ^ 

leading  to  the  two  separate  Helmholtz  equations: 


/(*)«  = -A/(jc).  (8) 

g{t)t,  =  - c2lg{t ).  (9) 


After  applying  boundary  conditions 

v2/M  =  -■*/(*)  l,,D- 


(10) 
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f(x)D, 


sin(kn  x), 


,  _ .  2  n  _  nn 

hi  ~  L  (11) 


where  A  represents  the  eigenvalue  and  A„  is  the  wave¬ 
length.  Similarly  for  time: 

^2 8{t)  =  ~(r  Xg(t)\3D,  C=J‘  (12) 

g{t)Dfl  =  \pt  sin(°3„  t),  (On  =  c  kn  (13) 

This  formulation  leads  to  the  well-known  ’’Initial  Value 
Problem  (IVP)”  from  standard  differential  equation  books. 
Namely,  if  a  function  is  known  completely  across  the  x- 
axis  at  t  —  0  =  tQ  and  its  time-derivative  is  also  known  for 
all  x  at  t0,  the  complete  Fourier  spectra  can  be  determined 
for  all  time.  To  illustrate  this  solution,  consider  the  case 
when  there  are  mixed  boundaries  in  both  space  and  time. 
Under  those  circumstances,  the  general  solution  is  formed 
by  allowing  both  Dirichlet  and  Neumann  solutions  to  ex¬ 
ist.  The  basis  functions,  f(x)  and  git  )  are  known  from  the 
underlying  differential  equations.  In  this  case,  the  general 
solution  takes  the  form  (neglecting  any  source  terms): 


Fig.  4  Dual  Time  Problem:  at  two  times,  t\  and  tz,  the  function 
F{x,t)  is  completely  known  along  the  x-axis. 


Pm  =  £  [ACng(t0)DJ(x)DJf(x)NmdQ  + 

n= 0  J 

+  f  [BCng(t0)NJ(x)Djl]f(x)NmdQ  + 

+  f  [AD»g(t0)DJ(x)Nii\f(x)NwdQ  + 
+  f  [BDng{tQ)Nnf{x)Nn}f{x)NmdQ 


F{x,t) 

Hx,t0) 


CCm 


X  [A»^(0An+B»5(?)jV,„]  [Cnf{.x)Dn+Dnf(x)Nl] 

n= 0 

X  [AC^(f0  )DJ{X)D,n  +BCng{t0)NJ{x)Dn]  + 

n= 0 

+  [AL>ng(t0)DJ(x)Nn  +  BDng  (t0  )Nn  f  (x)Nn  \ 

/  F{x,t0)dalJ(x)DmdQ 

X  /  [■ ACng(.t0)D,nf(x)D,n]f(x)D,»,d£2  + 
n=0  ■’ 

+  f  [BCng{t0)NJ{x)Dn]f{x)DmdQ  + 

+  f  [AF>ng{t0)DJ{x)Njl]f{x)DmdQ  + 

+  [ BDng(t0)Nnf(x)Nn]f(x)DmdQ 

X  [ACng(t0)D,n+BCng(t0)N,n\ 

n=0  J 

+  [AB>ng{t0)Dn  +BDng(t0)Nn]  j fNnfDmdQ 

X  \_A^ng^t0^D,n  + 

n= 0 


Pm  =  X  [' kCng(t0)Dn+BCng(t0)Nn \  f  fnJNmdQ 

n=0  ■’ 

+  [ADng(t0)Dn  +BDng(t0)Nn\  j)  fNnfNmdQ. 

Pm  =  X  [ACng<d 0)d,„  +BCng{t0)N,n\  0  + 

n— 0 

+  [ADng(t0)Dn  BDng(t0)Nnj  8nm 
P„  =  [AF>ng{t0)Dn  +BDng(t0)Nin] 

Pn  =  j  F(x,t0)f(x)NndQ 

The  problem  with  using  mixed  boundaries  becomes  ap¬ 
parent  quickly.  There  are  four  unknowns,  ACn.  BCn,ADn, 
and  BD„,  yet  only  two  equations  of  constraint  coming  from 
the  initial  conditions,  namely  from  F(x,t0),  we  get  an  and 
j3„.  Two  more  equations  are  needed  in  order  to  fully  de¬ 
termine  the  amplitudes  of  the  basis  functions.  The  Ini¬ 
tial  Value  Problem  uses  knowledge  of  two  pieces  of  in¬ 
formation,  F(x,t0)  and  ^F(x,t |?0).  To  set  up  the  problem 
attempted  on  the  Chesapeake  Bay,  the  IVP  is  altered  to 
a  new  form,  namely,  the  “Dual  Time  Problem”,  (DTP), 
where  a  second  time  is  introduced  instead  of  using  |f0). 


+  [AF>ng(,t0)Dn  +BD„g(t0)Nn]  0 
a»  =  [ACng{to)D,n  +BCng{t0)Nfl\ 
a„  =  j  F(xj0)  fix)  DndQ 
similarly  for  the  Neumann  term: 

Pm=  j F(x,t0)datJ(x)NmdQ 


4  Dual  Time  Problem 

Given  complete  knowledge  of  a  function  in  x  at  two  times, 
F (x,  tl )  and  Fix,  t2),  four  quantities  are  defined,  an ,  /3„ ,  yn , 
and  An.  From  these  four,  the  spectral  components  are  fully 
determined. 
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Fig.  5  Once  the  coefficients  have  been  calculated,  the  solution  may 
be  moved  forward  in  time,  allowing  for  prediction. 


««=  jF{x,tl)datJ{x)DndQ 
Pn  =  jF{x,tl)datJ{x)NndQ 
Yn=  j>  F(x,  t2)dala  f  (x)Dll dQ 
A„=  j)  F(x,t2  )data  f  {x)Nn  dQ 

an  —  \ACng{tl)Dn-\-BCng{t])Nn\ 

Pn  =  \ADn8{ti  )D,n  +  BDng(tl  )Ar  „] 

In  =  [- ACng(t2)Dn+BCng(t2)Nn ] 

A»  =  [- ADng(t2)Dn+BDng(t2)Ntn ] 

(  &n  Yn  \  (  ACn  BC„  \  (  g(t\)j)n  g(t 2)o,n\ 

V  Pn  An)  ~  \ADn  BDn)  \  g(U)Nn  g(t2)Nn  ) 

ACn  BC„  \  _  (  an  y,\  f  g(t i)D„  g(t2)D^  V1 
A D n  BDn  )  y  Pn  An  )  y  g(t\  )Nn  g{t2)NrI  ) 

Several  important  things  to  note  about  this  result.  This  is 
really  just  the  IVP  masked  over  two  times.  The  stablity 
of  the  solution  depends  on  having  good  data  covering  the 
whole  x-domain.  Because  the  last  matrix  in  the  derivation 
must  be  inverted,  it  cannot  be  singular  or  close  to  singu¬ 
lar.  As  a  result,  it  is  important  to  pick  two  times  t\  and 
ti  that  are  not  on  any  nodes  for  the  temporal  basis  func¬ 
tions.  Provided  these  conditions  are  met,  the  coefficients 
obtained  should  allow  the  function  F(x,t)  to  be  projected 
ahead  (or  behind)  in  time.  Finally,  because  the  basis  set 
was  chosen  in  time  to  be  sines  and  cosines,  there  is  tem¬ 
poral  repetition  on  the  time  scale  of  T  —  L  /  c.  Outside 
of  this  time  window,  the  temporal  components  repeat,  pre¬ 
venting  prediction  further. 


5  Dual  Position  Problem  (conjugate  to  time  problem) 

The  DTP  is  not  new,  nor  is  it  controversial;  it  is  simply  an 
extension  of  the  IVP.  The  goal  of  this  paper  is  to  present 


Fig.  6  Dual  Position  Problem:  at  two  locations,  xi  and  X2,  the  func¬ 
tion  F(x,t)  is  completely  known  along  the  f-axis. 


the  Dual  Position  Problem  (DPP).  Although  a  simple  trans¬ 
formation  from  the  time  axis  to  the  spatial  axis,  the  re¬ 
sults  are  contentious.  Simply  put,  by  collecting  data  at 
two  spatial  locations  over  a  long  enough  time  window. 
At  >=  Ax max  /  c ,  the  spectral  components  are  determined 
and  applied  to  the  spatial  domain,  allowing  spatial  predic¬ 
tion. 
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Fig.  7  Once  the  coefficients  have  been  calculated,  the  solution  may 
be  used  to  solve  across  the  entire  spatial  domain. 


ccn  =  ACJ(x\  )Dn  +ADJ{x\  )Nn 
Pn  —  BCJ{xi)Dn  +BDJ(x\)Nn 
Yn=ACJ{x1)Dn+ADJ(x1)Nn 
An  =  BCnf(x2)r,n+BDJ(x2)N,n 

f  (/-a  Yn  \  _  f  AC„  ADn  \  f  f(x  i)Dn 

\  Pn  An)  -  \BCnBDn)  \f{xx)Ntn  f{x2)NJ 

f  ACn  ADn  \  /  OC),  Yn  \  (  )  D.n 

\BCn  BDn)  -  ypn  An)  \f{xx)Nn  f(x2)Nn  ) 

Similar  to  the  DTP,  several  important  things  to  note 
about  this  result.  This  is  really  just  the  DTP  in  space,  not 
time.  The  stablity  of  the  solution  depends  on  having  good 
data  covering  the  whole  /-domain.  Because  the  last  matrix 
in  the  derivation  must  be  inverted,  it  cannot  be  singular  or 
close  to  singular.  As  a  result,  it  is  important  to  pick  two 
locations  x\  and  x2  that  are  not  on  any  nodes  for  the  spa¬ 
tial  basis  functions.  Provided  these  conditions  are  met,  the 
coefficients  obtained  should  allow  the  function  F(x,t)  to 
be  projected  in  space.  Finally,  because  the  spatial  basis 
set  was  computed  in  COMSOL,  confidence  in  the  solu¬ 
tion  is  paramount.  The  whole  derivation  is  hinged  on  or¬ 
thonormality,  so  when  COMSOL  computes  a  “tight”  set 
of  functions,  this  method  should  work.  For  most  complex 
boundaries,  it  is  not  practical  to  extend  the  solution  be¬ 
yond  the  stated  spatial  boundaries,  however,  by  choosing 
a  time  window  long  enough  to  span  the  largest  extant  of 
the  spatial  domain,  this  issue  should  not  be  a  problem. 


6  Results  and  Conclusions 

The  coefficients  computed  in  spectral  domain  are  assumed 
to  be  constant  in  both  space  and  time  for  the  windows  un¬ 
der  consideration.  Nature,  however,  is  not  static.  By  con¬ 
tinuing  to  take  data  at  spatial  locations  over  time,  overlap¬ 
ping  time  windows  may  be  employed  to  re-calculate  the 
coefficients.  In  the  new  time  window,  the  coefficients  will 
still  be  constants,  however,  they  will  represent  the  aver¬ 
age  spectral  coefficients  for  that  specific  timing  window. 


AC(n) 


Fig.  8  By  calculating  overlapping  windows  in  time,  the  spectra  cal¬ 
culated  are  shown  to  exhibit  time-dependence. 


By  successively  overlapping  time  windows  and  comput¬ 
ing  spectral  components,  a  time  series  for  the  spectra  is 
obtained,  as  illustrated  in  figure  8.  From  this  process,  time 
prediction  may  be  achieved  by  fitting  the  spectral  time 
variation  to  traditional  functions,  polynomials,  sinusoid, 
exponential,  etc... 

Working  with  a  toy  model  shed  some  light  on  this 
problem.  The  goal  was  to  take  a  user-defined  function  at 
two  times  along  the  x-axis,  use  the  DTP  to  calculated  the 
coefficients  so  that  the  solution  could  be  calculated  along 
the  time  axis.  Once  solved  along  the  time  axis,  the  test 
would  be  to  extract  the  solution  along  two  spatial  loca¬ 
tions,  x\  and  A'2,  and  re-calculate  the  spectra  using  the  DPP. 
If  the  solution  is  robust,  the  original  user-defined  functions 
should  be  reproducable  using  the  DPP  spectra  to  recon¬ 
struct  along  the  spatial  domain.  Early  attempts  to  perform 
this  robust  test  of  these  methods  failed. 

Careful  inspection  of  the  reconstructed  results  revealed 
that  the  methods  begin  to  falter  as  the  solution  is  pro¬ 
jected  towards  the  edges  of  the  domains,  both  time  and 
space.  The  reason  stems  from  the  fact  that  the  correct  so¬ 
lution  requires  a  summation  over  all  modes,  from  n  —  1 
up  to  n  —  <x>.  When  fewer  modes  are  used,  taking  Nmax 
from  the  Nyquist  limit,  the  solution  begins  to  fail  in  its  re¬ 
construction  as  it  approaches  the  far  edges  of  the  domain. 
Further  preventing  the  solution  from  matching  predictions 
beyond  the  windows  is  the  fact  that  harmonic  functions 
repeat  outside  of  thier  stated  domains,  so  a  fucntion  that 
would  normally  increase  linearly  in  time,  becomes  a  saw¬ 
tooth  function,  suddenly  dropping  as  t  approaches  its  lim¬ 
its.  This  abrupt  change  in  the  solution  prevents  the  toy 
model  from  easily  tranferring  spectra  from  the  time  do¬ 
main  to  the  spatial  domain  and  then  accurately  transform¬ 
ing  back  to  the  time  domain,  for  validation.  To  address  this 
issue  about  crossing  the  domain  windows,  a  linear  term  is 
added  to  both  the  time  and  space  functions.  The  end  re¬ 
sult  is  16  terms,  requiring  8  spatial  locations  to  fully  solve 
for  the  spectral  coefficients.  In  order  to  avoid  a  singular 
matrix  inversion,  by  choosing  10  spatial  locations  instead 
of  8,  when  a  particular  matrix  is  close  to  being  singular, 
the  spectra  obtained  tend  to  vary  significantly.  Rather  than 
throw  out  these  results,  for  the  set  of  10  locations,  every 
set  of  8  locations  is  used  to  calculate  the  spectra,  giving 
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Fig.  9  Realistically,  more  than  two  positions  will  be  used  to  fully 
calculate  the  coefficients,  leading  to  better  prediction  in  both  time 
and  space. 


results  which  are  pushing  our  traditional  approaches  to 
older  problems,  such  as  signal  processing.  Problems  en¬ 
countered  with  the  Chesapeake  have  simply  been  unprece¬ 
dented  in  the  past  due  to  a  lack  of  confidence  in  the  so¬ 
lutions  computed  for  the  spatial  domain.  Through  COM¬ 
SOL,  one  must  re-think  the  entire  avenue  set  down  by  the 
one  dimensional  (time-based)  signal  processing  commu¬ 
nity. 
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a  set  of  results  for  the  coefficients  of  size  \N\/S\  —  45, 
where  N  is  the  total  number  of  locations  taken.  Taking  the 
mean  of  the  spectral  results  prevents  any  one  result  from 
skewing  the  data  set. 
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The  Chesapeake  Bay  remains,  and  will  so  for  a  while, 
a  work  in  progress.  This  latest  revelation,  that  one  can 
use  a  small  number  of  locations  to  monitor  the  entire  Bay 
by  taking  sufficient  time  windows  presents  a  hopeful  out¬ 
look  for  future  signal  characterization  of  the  Bay.  The  util¬ 
ity  to  the  COMSOL  community  is  more  towards  appli¬ 
cation.  By  using  a  high-performance  computing  applica¬ 
tion  like  COMSOL,  scientists  are  capable  of  producing 
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